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LONG-TERM  GOALS 

Our  goal  is  to  create  a  robust  and  highly  efficient  nonhydrostatic  code  that  can  be  applied  in  a 
predictive  manner  to  littoral  zones  and  bays  on  a  spatial  scale  of  the  order  of  100  km  and  time  scales  of 
the  order  of  days  to  weeks.  The  product  will  be  a  user-oriented  production  code  for  the  coastal 
environment. 

OBJECTIVES 

This  project  aims  to  synthesize  a  number  of  computational  and  numerical  tools  to  produce  an 
innovative  and  powerful  coastal  ocean  simulation  tool  to  generate  accurate  predictions  of  motions  and 
transport  in  coastal  oceans  under  conditions  when  a  nonhydrostatic  and  terrain-following 
representation  is  essential.  Another  primary  objective  is  to  apply  the  best  algorithms  available  to  create 
a  simulation  code  that  is  capable  of  representing  the  physical  processes  at  high  resolution  and  capable 
of  very  high  speed  on  multiple-processor  parallel  computers. 

APPROACH 

This  project  is  complementary  to  NSF  ITR  Grant  OCE-01 131 1 1.  That  project  is  focused  on  creation  of 
a  new  generation  computer  code  and  application  of  it  to  internal  waves  in  Monterey  and  Mamala  Bays. 
Four  faculty  are  involved.  Robert  Street  is  the  PI.  Margot  Gerritsen  is  an  expert  on  numerical  methods 
and  adaptive  mesh  refinement  strategies.  Mark  Merrifield  is  a  member  of  the  HOME  project  and  will 
provide  the  data  for  Mamala  Bay  and  assist  in  its  interpretation.  The  main  code  development  is 
supported  by  the  NSF  ITR  project  and  Oliver  Fringer  is  carrying  out  the  development  work.  This  ONR 
funding  is  supporting  a  doctoral  student  to  take  responsibility  for  applications.  This  entails  the  crucial 
task  of  generation  of  grids  and  initial  data  and  boundary  forcing.  He  is  carrying  out  and  analyzing 
production  runs  as  well. 

We  will  employ  numerical  large-eddy  simulation  to  generate  accurate  predictions  of  motions  and 
transport  in  coastal  oceans  under  conditions  when  a  nonhydro static  and  terrain-following 
representation  is  essential.  The  simulation  should  be  able  to  handle  all  processes  from  the  free-surface 
to  the  bed  so  that  the  whole  range  from  surface  waves  through  internal  waves  to  sediment  motions  can 
be  simulated.  Two  specific  applications  of  the  completed  code  are  planned.  Both  involve 
nonhydrostatic  evolution  of  internal  tides.  The  first  is  in  Monterey  Bay,  California,  where  solitons  are 
expected  to  form  and  the  internal  tide  signal  is  intensified  in  the  bottom  of  the  Monterey  Canyon.  The 
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second  is  in  Mamala  Bay,  Hawaii,  where  high-amplitude  and  nonlinear  internal  tides  are  observed  and 
there  are  key  questions  to  be  answered  about  the  nature  of  these  waves  and  their  behavior  in  the  Bay. 

From  a  numerical  point  of  view  our  approaches  include: 

1 .  Parallel  processing:  We  are  using  MPI,  the  message  passing  interface. 

2.  Free-surface  and  nonhydrostatic  Poisson-equation  solvers:  We  are  exploring  the  use  of  PETSc 
[Portable  Extensible  Toolkit  for  Scientific  Computing ;  see  below]. 

3.  Accurate  advection:  We  are  employing  mass,  volume,  momentum  and  energy  conserving 
Eulerian  schemes. 

4.  Large  Eddy  Simulation:  We  have  developed  and  tested  a  mixed  subfdter-scale  model  for 
turbulence. 

5.  State  of  the  art  gridding:  Produces  high  quality  unstructured  grids  in  the  horizontal  and 
structured  in  the  vertical  for  optimal  compromise  between  efficiency  and  flexibility/grid  quality. 

WORK  COMPLETED 

For  Fiscal  Year  2003,  the  primary  tasks  under  this  ONR  grant  were  the  analysis  and  preparation  of 
highly  resolved  unstructured  grids  for  Monterey  Bay  and  to  begin  production  runs  of  the  parallel  code 
to  study  internal  waves  in  idealized  shelf  geometries  and  Monterey  Bay.  Currently,  we  are  using  a 
three-dimensional  primitive  equation  version  of  our  unstructured-grid  code  and  are  testing  its  ability  to 
(1)  generate  internal  waves  at  the  shelf  break  in  simplified  geometries  and  (2)  compute  the  hydrostatic 
internal  wave  spectrum  in  Monterey  Bay.  The  code  employs  the  momentum  and  energy  conserving 
Eulerian  formulation  of  Perot  (2000)  for  advection  of  momentum,  and  the  mass  conserving 
formulation  of  Gross  et  al.  (2002)  for  scalar  advection.  The  Poisson  equations  resulting  from  the  semi- 
implicit  free-surface  fonnulation  (Casulli,  1998)  and  the  nonhydrostatic  pressure  are  being  solved  with 
the  preconditioned  conjugate  gradient  algorithm.  Prof.  Gerritsen  is  exploring  the  use  of  PETSc  to 
parallelize  and  optimize  these  solvers. 

RESULTS 

Advection  schemes 

A  significant  amount  of  work  this  year  has  focused  on  the  choice  of  the  advection  scheme  for 
momentum.  Currently  we  are  employing  a  momentum  and  energy  conserving  Eulerian  formulation 
(Perot,  2000)  for  unstructured,  staggered  grids.  The  original  intent  was  to  employ  the  Euler-Lagrange 
formulation  of  Casulli  (1998)  with  Kriging  (LeRoux  et  al.  1997)  as  the  interpolator.  However,  because 
the  parallel  fonnulation  limits  the  Courant  number  to  less  than  unity  to  prevent  Euler-Lagrange 
tracebacks  from  moving  into  a  neighboring  processor,  it  was  decided  that  an  Eulerian  formulation 
would  be  more  appropriate. 

Treatment  of  the  bottom  boundary 

We  have  also  analyzed  the  treatment  of  the  bottom  boundary  carefully  in  order  to  determine  an 
appropriate  means  of  accurately  representing  the  bottom  bathymetry.  The  current  fonnulation  employs 
the  semi-implicit  free-surface  formulation  of  Casulli  (1998),  which  requires  a  stair-stepped  treatment 
of  the  bottom  bathymetry.  While  this  formulation  removes  the  false  baroclinic  pressure  gradients 
associated  with  sigma  coordinate  grids,  the  trade  off  introduces  errors  in  along-slope  gradients  which 
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are  not  as  well  resolved  with  stair-stepped  grids.  The  shaved  cell  approach  of  Adcroft  et  al.  (1997)  was 
considered,  but  this  approach,  along  with  other  similar  cut-cell  techniques  (e.g.  Kirkpatrick  et  ah, 

2003)  is  not  feasible  when  wetting  and  drying  is  desired.  The  original  intent  was  to  employ  the 
immersed  boundary  method  (IBM),  which  has  been  employed  by  Tseng  and  Ferziger  (20003)  for 
geophysical  applications  and  showed  improvement  over  stair-stepped  grids.  Despite  their  favorable 
results,  we  have  continued  to  use  our  stair-stepped  formulation  because,  as  opposed  to  the  IBM 
method,  it  guarantees  conservation  of  mass  and  volume,  and  we  are  wary  of  the  use  of  non¬ 
conservative  schemes  for  the  high  Reynolds  number  flows  encountered  in  our  application  where  a  drag 
law  must  be  used  to  model  the  shear  stress  at  the  bottom  bathymetry.  We  will  continue  to  search  for 
improved  conservative  bottom-following  formulations. 

Grid  generation 


Unstructured  grid  generation  has  proven  to  be  a  complex  task  because  a  triangulation  in  which  the 
elements  are  generated  to  follow  or  resolve  certain  regions  of  the  complex  bathymetry  of  Monterey 
Bay  always  results  in  degenerate  triangles.  Triangles  are  considered  degenerate  when  the  Voronoi 
points  (which  are  the  centers  of  the  circumcircles  of  the  triangle  vertices)  lie  outside  of  the  triangles 
because  of  the  creation  of  obtuse  triangles.  While  it  is  possible  to  minimize  the  frequency  of 
degenerate  triangles  by  specifying  the  minimum  angle  allowed  in  a  triangulation,  the  Delaunay 
triangulation  of  Shewchuck  (2002)  is  guaranteed  to  converge  only  when  the  minimum  angle  is  less 
than  33.8°.  As  an  example,  Figure  1  depicts  a  triangulation  of  Monterey  bay  in  which  triangle  faces 
are  coincident  with  chosen  isobaths.  From  the  figure,  it  is  clear  that  an  arbitrary  choice  of  isobaths 
results  in  degenerate  triangles,  especially  when  isobaths  are  very  close  together,  in  which  case  large 
area  ratios  between  neighboring  triangles  are  also  generated. 


Figure  1.  Isobaths  at  0,  345, 1122, 1 760,  2416,  and  2796  m  and  the  associated  Voronoi 
triangulation  for  which  the  isobaths  coincide  with  triangle  faces. 

An  alternative  is  to  employ  BatTri  (Bilgili  and  Smith,  2001),  a  front-end  for  Triangle,  which  allows  the 
user  to  select  regions  of  the  bathymetry  for  improved  refinement.  Examples  of  using  BatTri  are  shown 
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in  Figure  2.  It  can  be  seen  that,  while  refinement  is  achieved  in  a  straight-forward  manner,  the 
transition  of  grid  area  across  the  refinement  boundary  is  much  too  high.  While  Triangle  does  allow  the 
minimum  angle  constraint,  it  does  not  allow  a  constraint  on  the  maximum  area  ratio.  As  it  is,  Triangle 
can  be  used  as  a  library  function  directly  in  the  main  code,  but  the  BatTri  routines  cannot  be  employed 
directly  because  they  are  a  front-end  to  Triangle  and  are  written  in  Matlab.  We  are  considering  the 
GAMBIT  (Fluent,  Inc,  2001)  unstructured  grid  generating  package,  which  will  hopefully  resolve  these 
issues. 


Figure  2.  Examples  of  using  the  front  end  to  BATTRI,  showing  the  result  of  refining  a  region  of  the 
bay  as  well  as  imposing  a  constraint  on  the  triangulation  to  follow  the  canyon  bathymetry 


1  meter  M2  tide  (period  of  12.42  hrs) 


100  km 


500  m 


Figure  3.  The  simplified  two-dimensional  shelf-break  bathymetry,  showing  the  forcing  and 
buoyancy  frequencies.  The  effect  of  varying  the  slope  angle  a  is  shown  in  Figure  4. 

Internal  waves  in  idealized  geometries 

To  verify  the  capability  of  the  code  to  generate  internal  waves  at  a  shelf  break,  we  have  tested  it  using 
the  simplified  coastal  shelf  bathymetry  shown  in  Figure  3.  The  stratification  is  linear  with  a  constant 
buoyancy  period  27t/N  of  3  hours.  The  results  of  forcing  the  left  boundary  with  an  M2  tidal  amplitude 
of  1  m  (2n/(»=  1 2.42  hours)  and  four  different  slope  angles  a  are  shown  in  Figure  4.  Consistent  with 
the  linearized  dispersion  relation  which  determines  the  propagation  angle  of  the  group  velocity  vector, 
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the  internal  waves  propagate  away  from  the  shelf  break  at  an  angle  given  by  sin'1(co/N)=14.3°.  The 
four  different  shelf-break  slopes  in  Figure  4  are  given  by  9.5°  (subcritical),  14.3°  (critical),  26.6° 
(supercritical),  and  90°  (vertical).  The  seaward  propagating  depth-integrated  energy  flux  is  shown  for 
each  case  as  a  multiple  of  the  energy  flux  for  the  critical  case,  Ec.  The  results  in  Figure  4  show  that  the 
energy  flux  increases  as  the  slope  angle  increases  beyond  the  critical  slope.  This  results  from  the 
increase  in  the  vertical  velocity  at  the  shelf  break  as  the  slope  angle  increases,  which  is  directly  related 
to  the  strength  of  the  internal  wave  generation  zone  at  the  shelf  break.  Internal  wave  generation  is 
weakest  for  the  subcritical  slope,  but  energy  flux  is  also  reduced  because  the  subcritical  slope  impedes 
seaward  propagation  of  internal  wave  energy.  Numerical  oscillations  are  pronounced  in  the  vertical 
slope  case  because  of  the  sharp  discontinuity  in  the  bathymetry. 


Subcritical 


E 


Critical 


J 


1-3Ec  ◄ 


Figure  4.  Perturbation  density  fields  resulting  from  the  generation  of  an  internal  tide  at  a  shelf 
break,  showing  the  depth-integrated  energy  flux  for  each  case. 


IMPACT/APPLICATIONS 


A  key  goal  is  to  produce  a  production-style  code,  i.e.,  one  that  can  be  used  by  others,  is  fully 
documented  and  tested,  and  in  the  public  domain.  Thus,  a  result  will  be  a  robust  tool  for  use  by  the 
coastal  oceanographic  community  to  study  coastal  processes. 

Previous  studies  have  documented  the  importance  of  gaining  an  understanding  of  the  generation, 
propagation,  and  fate  of  internal  tidal  waves  on  the  continental  shelf  and  coastal  bays.  Definitive 
answers  about  the  role  of  these  waves  would  materially  increase  our  understanding  of  the  mixing  and 
transport  on  the  shelves  and  in  the  bays  and  their  influence  on  sediment  transport.  Flence  a  predictive 
capability  has  broad  coastal  oceanographic  applications,  including  for  example,  to  the  modeling  phase 
of  the  HOME  project. 
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RELATED  PROJECTS 


NSF  OCE-Ol  13 1 1 1  High  Resolution  Coastal  Ocean  Modeling.  The  numerical  code  for  this  ONR 
project  is  being  created  under  the  NSF  grant. 
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